#Import rcas and reaches

fgdb <- "Spatial_data/Anchor/Anchor.gdb"

# Read the rca and reaches feature classes
rcas <- sf::st_read(dsn = fgdb, layer = "anchor_rcas_attributed_06022020") %>% select(c(rca_id, Shape, Shape_Area,rca_elev_mn))
## Reading layer `anchor_rcas_attributed_06022020' from data source `W:\GitHub\Temperature_Data\Spatial_data\Anchor\Anchor.gdb' using driver `OpenFileGDB'
## Simple feature collection with 1210 features and 25 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 119365 ymin: 1077505 xmax: 158294.3 ymax: 1108700
## epsg (SRID):    3338
## proj4string:    +proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
rca_reach <- sf::st_read(dsn = fgdb, layer = "anch_rca_reaches_attributed_06022020") %>% mutate(rca_id=reach_id) %>% select(c(rca_id, Shape, Shape_Length))
## Reading layer `anch_rca_reaches_attributed_06022020' from data source `W:\GitHub\Temperature_Data\Spatial_data\Anchor\Anchor.gdb' using driver `OpenFileGDB'
## Simple feature collection with 1212 features and 8 fields
## geometry type:  MULTILINESTRING
## dimension:      XY
## bbox:           xmin: 119471.1 ymin: 1077726 xmax: 157727.9 ymax: 1108389
## epsg (SRID):    3338
## proj4string:    +proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

#load chinook predictions dataframe and join to spatial data

#load chinook predictions dataframes
load(file = "output/K_p.Rdata")
K_p2 <- K_p %>% select(rca_id, K_p, K_s, K_r, preds, year, month, day, date, decade)
                       
#helpful for how to join dataframe to spatial data maintaining sf class: #http://www.nickeubank.com/wp-content/uploads/2015/10/RGIS2_MergingSpatialData_part1_Joins.html

K_p_rca <- sp::merge (rcas, K_p2, by="rca_id") 

#remove large objects from global environment
remove(K_p2, K_p)

#Create dataframes for maps and delete those that are really large objects

#filter for one day to test map template
K_try <- K_p_rca %>% filter(day==202 & year==2004)

#create datasets needed for plot -  interannual variability of rcas
annual <- K_p_rca %>% group_by(rca_id, year) %>% summarize(yr_mn=mean(preds), yr_sd=sd(preds))

#faceted map for cold year and warm year - annual average temperatures
warm_cold <- K_p_rca %>% filter(year==2004 | year==1992) %>% group_by(rca_id, year) %>% 
  summarize(yr_mn=mean(preds), yr_sd=sd(preds))

#average temperatures by decade and mapped
decade<- K_p_rca %>%  group_by(rca_id, decade) %>% summarize(dec_mn=mean(preds), dec_sd=sd(preds))

#average temperatures by decade and month
mth_decade<- K_p_rca %>%  group_by(rca_id, month, decade) %>% 
  summarize(dec_mn=mean(preds), dec_sd=sd(preds))

#animated plot1
K_aug <- K_p_rca %>% filter(month==8) %>% group_by (rca_id, year, month) %>% summarise(aug_mn = mean(preds))

#animation for whole year
K_2004 <- K_p_rca  %>% filter(year==2004) %>% group_by (rca_id, date) %>% summarise(day_mn = mean(preds)) 

#remove large objects from global environment
remove(K_p_rca)

#Static plots

#let's try a basic map of mean elevation of rcas with color gradient
ggplot(data=rcas) + geom_sf(aes(fill=rca_elev_mn)) + scale_colour_viridis_c(option = "C")

#filter predictions for one day and try creating base map of predictions with color gradient and change breaks in color ramp

ggplot(data=K_try) + geom_sf(aes(fill=preds)) + 
    scale_fill_gradient(low = "blue", high = "red", breaks=c(4,5,6,7,8,9,10,11,12,13,14,15))

#Map of inter-annual variability of predictions
ggplot(data=annual) + geom_sf(aes(fill=yr_sd)) +
    scale_fill_gradient(low = "blue", high = "red",name = "SD of Temperature (Celsius)") +
    labs(title="Interannual variation of mean annual temperatures", x="Longitude", y="Latitude") 

ggsave("Anchor_chinook_plots/Interannual_variation_mean_annual_temperatures.jpeg", plot=last_plot())
## Saving 7 x 5 in image
#warm and cold year compared
ggplot(data=warm_cold) + geom_sf(aes(fill=yr_mn)) + 
    scale_fill_gradient(name = "Temperature (Celsius)", low = "blue", high = "red") + facet_wrap(~year) +
     labs(title="Mean Annual Temperatures for warmest and coldest years", x="Longitude", y="Latitude")

ggsave("Anchor_chinook_plots/Mean_Annual_Temperature_warm_coldyears.jpeg", plot=last_plot())
## Saving 7 x 5 in image
#Mean decadal temperatures
ggplot(data=decade) + geom_sf(aes(fill=dec_mn)) + 
    scale_fill_gradient(name = "Temperature (Celsius)", low = "blue", high = "red") + facet_wrap(~decade) +
   labs(title="Mean decadal temperatures", x="Longitude", y="Latitude")

ggsave("Anchor_chinook_plots/Mean_decadal_temperatures.jpeg", plot=last_plot())
## Saving 7 x 5 in image

#Monthly averages faceted by decade - static

jun <- mth_decade %>% filter(month==6)
ggplot(data=jun) + geom_sf(aes(fill=dec_mn)) + 
    scale_fill_gradient(name = "Temperature (Celsius)", low = "blue", high = "red") + facet_wrap(~decade) +
   labs(title="Mean June temperatures by decade", x="Longitude", y="Latitude")

ggsave("Anchor_chinook_plots/Mean_june_decades.jpeg", plot=last_plot())
## Saving 7 x 5 in image
jul <- mth_decade %>% filter(month==7)
ggplot(data=jul) + geom_sf(aes(fill=dec_mn)) + 
    scale_fill_gradient(name = "Temperature (Celsius)", low = "blue", high = "red") + facet_wrap(~decade) +
   labs(title="Mean July temperatures by decade", x="Longitude", y="Latitude")

ggsave("Anchor_chinook_plots/Mean_july_decades.jpeg", plot=last_plot())
## Saving 7 x 5 in image
aug <- mth_decade %>% filter(month==8)
ggplot(data=aug) + geom_sf(aes(fill=dec_mn)) + 
    scale_fill_gradient(name = "Temperature (Celsius)", low = "blue", high = "red") + facet_wrap(~decade) +
   labs(title="Mean August temperatures by decade", x="Longitude", y="Latitude")

ggsave("Anchor_chinook_plots/Mean_aug_decades.jpeg", plot=last_plot())
## Saving 7 x 5 in image
sep <-  mth_decade %>% filter(month==9)
ggplot(data=sep) + geom_sf(aes(fill=dec_mn)) + 
    scale_fill_gradient(name = "Temperature (Celsius)", low = "blue", high = "red") + facet_wrap(~decade) +
   labs(title="Mean September temperatures by decade", x="Longitude", y="Latitude")

ggsave("Anchor_chinook_plots/Mean_sep_decades.jpeg", plot=last_plot())
## Saving 7 x 5 in image

#Animations

#Animate timeseries maps
#https://d4tagirl.com/2017/05/how-to-plot-animated-maps-with-gganimate


#Now let's try small time-series animation by first averaging over August and showing mean predictions over all years

#https://ditheringdata.netlify.com/2018/01/01/gganimate/
#https://semba-blog.netlify.com/10/29/2018/animating-oceanographic-data-in-r-with-ggplot2-and-gganimate/
#https://aberdeenstudygroup.github.io/studyGroup/lessons/Visualising%20Spatial%20Data/Visualising_Spatial_Data_in_R/

m1 <- ggplot() + geom_sf(data=K_aug, aes(fill=aug_mn)) + 
    scale_fill_gradient(name = "Temperature (Celsius)", low = "blue", high = "red") + 
    labs(title= "Average August Stream Temperature: {as.integer(frame_time)}", x="Longitude", y="Latitude") +
     transition_time(year)

animate(m1, nframes=38, fps=1)

#how about trying daily temperatures for 1 chosen year

m2 <- ggplot() + geom_sf(data=K_2004, aes(fill=day_mn)) + 
    scale_fill_gradient(name = "Temperature (Celsius)", low = "blue", high = "red") + 
    labs(title= "Average Daily Stream Temperature: {frame_time}", x="Longitude", y="Latitude") +
     transition_time(date)

animate(m2, fps=1)

anim_save(filename="Daily_2004_animation.mp4", animation = last_animation())

```